Reconciling dark energy models with f(R) theories 
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Higher order theories of gravity have recently attracted a lot of interest as alternative candidates to 
explain the observed cosmic acceleration without the need of introducing any scalar field. A critical 
ingredient is the choice of the function f(R) of the Ricci scalar curvature entering the gravity 
Lagrangian and determining the dynamics of the universe. We describe an efficient procedure 
to reconstruct f(R) from the Hubble parameter H depending on the redshift z. Using the metric 
formulation of f(R) theories, we derive a third order linear differential equation for f(R(z)) which can 
be numerically solved after setting the boundary conditions on the basis of physical considerations. 
Since H(z) can be reconstructed from the astrophysical data, the method we present makes it 
possible to determine, in principle, what is the f(R) theory which best reproduces the observed 
cosmological dynamics. Moreover, the method allows to reconcile dark energy models with f(R) 
theories finding out what is the expression of f(R) which leads to the same H(z) of the given 
quintessence model. As interesting examples, we consider "quiessence" (dark energy with constant 
equation of state) and the Chaplygin gas. 

PACS numbers: 98.80.-k, 04.50+h, 98.80.Jk, 98.80.Es 



I. INTRODUCTION 

The impressive amount of astrophysical data which 
have been accumulated in recent years has depicted a 
new standard cosmological model according to which the 
universe is spatially flat and undergoing a phase of ac- 
celerated expansion. Strong evidences in favour of this 
scenario are the Hubble diagram of Type la Supernovae 
(hereafter SNela) Q , the anisotropy spectrum of the cos- 
mic microwave background radiation (hereafter CMBR) 
and the matter power spectrum as measured by the 
clustering properties of the large scale distribution of 
galaxies [3j and the data coming from the Lya clouds 
01- Moreover, the abundance of clusters of galaxies and 
the gas mass fraction in clusters [j| constrain the matter 
density parameter Qm ~ 0.3 thus giving rise to the need 
for a new component with negative pressure to both close 
the universe (fltot — 1) and drive its accelerated expan- 
sion. This is what is usually referred to as dark energy 
whose subtle and elusive nature has opened the doors 
to an overwhelming flood of papers presenting a great 
variety of models which try to explain this phenomenon. 

The simplest explanation claims for the cosmological 
constant A thus leading to the so called ACDM model |6| ■ 
Although being the best fit to most of the available as- 
trophysical data , the ACDM model is also plagued by 
many problems on different scales. If interpreted as vac- 
uum energy, A is up to 120 orders of magnitudes smaller 
than the predicted value. Furthermore, one should also 
solve the coincidence problem, i.e. the nearly equivalence 
of the matter and A contribution to the total energy den- 
sity. As a response to these problems, much interest has 
been devoted to models with dynamical vacuum energy, 
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dubbed quintessence 0, . These models typically in- 
volve scalar fields with a particular class of potentials, 
allowing the vacuum energy to become dominant only 
recently (see E3 f° r comprehensive reviews). Although 
quintessence by a scalar field is the most studied candi- 
date for dark energy, it generally does not avoid ad hoc 
fine tuning to solve the coincidence problem. 

On the other hand, it is worth noting that, despite 
the broad interest in dark matter and dark energy, their 
physical properties are still poorly understood at a fun- 
damental level and, indeed, it has never been shown that 
they are actually two different ingredients apart the ob- 
servational fact that dark matter furnishes evidences of 
its existence in clustered structures (as galaxies, clusters 
of galaxies, etc.) and dark energy is an unclustered com- 
ponent which acts to accelerate the cosmic Hubble flow 
[llL Ha ] . These considerations motivated the great inter- 
est recently devoted to completely different approaches 
to the quintessence problem. Rather than fine tuning a 
scalar field potential (which, in several cases, is nothing 
else but a phenomenological ingredient not motivated by 
fundamental theories), it is also possible to explain the 
acceleration of the universe by introducing a cosmic fluid 
with an equation of state causing it to act like dark mat- 
ter at high densities and dark energy at low densities. 
Also this approach is, in some sense, phenomenological 
but it shows the attractive feature that these models can 
explain both dark energy and dark matter by a single 
mechanism (thus automatically solving the coincidence 
problem) and have therefore been referred to as unified 
dark energy (UDE) or unified dark matter (UDM). Some 
interesting examples are the generalized Chaplygin gas 
13 1, the tachyonic fiel d 111 , the condensate cosmology 
15 1 , the Hobbit model [16| , and the scaling dark energy 
17 1 . It is worth noting, however, that such models seems 
to be seriously affected by problems with structure for- 
mation 

El 

so that some work is still needed before they 
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can be considered as reliable alternatives to dark energy 
and dark matter. 

Actually, there is still a different way to face the prob- 
lem of cosmic acceleration. As stressed by Lue et al. 
[l9j . it is possible that the observed acceleration is not 
the manifestation of another ingredient in the cosmic pie, 
but rather the first signal of a breakdown of our under- 
standing of the laws of gravitation. From this point of 
view, it is thus tempting to modify the Friedmann equa- 
tions to see whether it is possible to fit the astrophysical 
data with a model comprising only the standard mat- 
ter. Interesting examples of this kind are the Cardassian 
expansion [23] and the DGP gravity pH ], 

In this same framework, there is also the attractive 
possibility to consider the Einstein General Relativity as 
a particular case of a more fundamental theory. This 
is the underlying philosophy of what are referred to as 
f(R) theories of gravity HI HI . Such theories 

are a class of the so called Extended Theories of Grav- 
ity which have become a sort of paradigm in the study 
of the gravitational interaction based on corrections and 
enlargements of the traditional Einstein - Hilbert scheme. 
In fact, in the last thirty years, several shortcomings 
came out in the Einstein scheme and people began to 
investigate whether General Relativity is the only fun- 
damental theory capable of explaining the gravitational 
interaction. Such issues come from cosmology and quan- 
tum field theory and essentially are due to the lack of 
a definitive quantum gravity theory. Alternative theo- 
ries of gravity have been pursued in order to attempt, 
at least, a semi - classical scheme where General Rela- 
tivity and its positive results could be recovered. The 
paradigm consists in adding higher - order curvature in- 
variants and minimally or non - minimally coupled scalar 
fields into dynamics which come out from the effective 
action of quantum gravity |27| . Other motivations to 
modify General Relativity come from the issue of a whole 
recovering of Mach principle (28) which leads to assume a 
varying gravitational coupling. All these approaches are 
not the "full quantum gravity but are needed as work- 
ing schemes toward it. In any case, they are going to 
furnish consistent and physically reliable results. Fur- 
thermore, every unification scheme as Superstrings, Su- 
pergravity or Grand Unified Theories, takes into account 
effective actions where non - minimal couplings or higher 
order terms in the curvature invariants come out. Such 
contributions are due to one - loop or higher - loop cor- 
rections in the high - curvature regimes near the full (not 
yet available) quantum gravity regime |27|. Specifically, 
this scheme is adopted in order to deal with the quanti- 
zation on curved spacetimes and the result is that the in- 
teractions among quantum scalar fields and background 
geometry or the gravitational self - interactions yield cor- 
rective terms in the Einstein - Hilbert Lagrangian [2^| . 
Moreover, it has been realized that such corrective terms 
are inescapable if we want to obtain the effective action 
of quantum gravity on scales closed to the Planck length 
|30| . Higher - order terms in curvature invariants (such 



as R 2 , R^R^, R^ aP R^ a p, ROR, or RD k R) or non- 
minimally coupled terms between scalar fields and geom- 
etry (such as <f) 2 R) have to be added to the effective La- 
grangian of gravitational field when quantum corrections 
are considered. Furthermore, from a conceptual point of 
view, there would be no a priori reason to restrict the 
gravitational Lagrangian to a linear function of the Ricci 
scalar R, minimally coupled with matter |3l| . Finally, 
the idea that there are no exact laws of physics but that 
the Lagrangians of physical interactions are stochastic 
functions - with the property that local gauge imparl- 
ances (i.e. conservation laws) are well approximated in 
the low energy limit and that physical constants can vary 
- has been taken into serious consideration . 

Besides fundamental physics motivations, all these the- 
ories have acquired a huge interest in cosmology due 
to the fact that they naturally exhibit inflationary be- 
haviours able to overcome the shortcomings of Standard 
Cosmological Model (based on General Relativity). The 
related cosmological models seem very realistic and, sev- 
eral times, capable of matching with the observations 
[33L l34| . Furthermore, it is possible to show that, via 
conformal transformations, the higher -order and non- 
minimally coupled terms always correspond to Einstein 
gravity plus one or more than one minimally coupled 
scalar fields [351 [IE [Sill]- More precisely, every higher - 
order term always appears as a contribution of order two 
in the equations of motion. For example, a term like 
R 2 gives fourth order equations [3j|, R OR gives sixth 
order equations [H,|4(]], RD 2 R gives eighth order equa- 
tions |41| and so on. By a conformal transformation, any 
2nd - order of derivation corresponds to a scalar field: for 
example, fourth - order gravity gives Einstein plus one 
scalar field, sixth order grav ity gives Einstein plus two 
scalar fields and so on |38l |42| . This feature results 
very interesting if we want to obtain multiple inflationary 
events since an early stage could select "very" large - scale 
structures (clusters of galaxies today), while a late stage 
could select "small" large-scale structures (galaxies to- 
day) 0] • The philosophy is that each inflationary era is 
connected with the dynamics of a scalar field. Further- 
more, these extended schemes naturally could solve the 
problem of " graceful exit" bypassing the shortcomings of 
former inflationary models |34l l43j . 

However, in the weak -field limit approximation, these 
theories are expected to reproduce General Relativity 
which, in any case, is experimentally tested only in this 
limit |44| . This fact is matter of debate since several rela- 
tivistic theories do not reproduce exactly Einstein results 
in the Newtonian approximation but, in some sense, gen- 
eralize them. As it was firstly noticed by Stelle |45|, a 
R 2 theory gives rise to Yukawa -like corrections to the 
Newtonian potential which could have interesting phys- 
ical consequences. For instance, some authors claim to 
explain the flat rotation curves of galaxies by using such 
terms . Others 0] have shown that a conformal the- 
ory of gravity is nothing else but a fourth - order theory 
containing such terms in the Newtonian limit. Besides, 
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indications of an apparent, anomalous, long range accel- 
eration revealed from the data analysis of Pioneer 10/11, 
Galileo, and Ulysses spacecrafts could be framed in a 
general theoretical scheme by taking corrections to the 
Newtonian potential into account |48j. In general, any 
relativistic theory of gravitation can yield corrections to 
the Newton potential (see for example which, in the 
post - Newtonian (PPN) formalism, could furnish tests 
for the same theory 01 • Furthermore the newborn grav- 
itational lensing astronomy |50j is giving rise to addi- 
tional tests of gravity over small, large, and very large 
scales which very soon will provide direct measurements 
for the variation of Newton coupling G n . the po- 
tential of galaxies, clusters of galaxies and several 
other features of gravitating systems. Such data will be, 
very likely, capable of confirming or ruling out the physi- 
cal consistency of General Relativity or of any Extended 
Theory. 

In this paper, we will restrict to f(R) theories and 
we will face the problem to find the cosmological mod- 
els (derived from a generic f(R)) consistent with data. 
In this case, the Fricdmann equations have to be given 
away in favour of a modified set of cosmological equations 
which are obtained by varying a generalized gravity La- 
grangian where the scalar curvature R has been replaced 
by a generic function f(R). The usual General Relativity 
is recovered in the limit f{R) = R, while completely dif- 
ferent results may be obtained for other choices of f(R). 
While in the weak field limit, in particular at Solar Sys- 
tem scales, the theory should give the usual Newtonian 
gravity, at cosmological scales there is an almost com- 
plete freedom in the choice of f(R) thus leaving open the 
way to a wide range of models 1 . The key point of f(R) 
cosmological models is the presence of modified Fried- 
mann equations obtained by varying the generalized La- 
grangian. However, the main problem of this approach 
lies here since it is not clear how the variation has to be 
performed. Actually, once the Robertson - Walker metric 
has been assumed, the equations governing the dynamics 
of the universe are different depending on whether one 
varies with respect to the metric only or with respect to 
the metric components and the connections. It is usual 
to refer to these two possibilities as the metric and the 
Palatini approach respectively. The two methods give 
the same result only in the case f(R) = R, while lead, in 
general, to significantly different field equations for other 
choices of the gravity Lagrangian (see [ME1II1E1 and 
references therein). 

Although much interest has been recently devoted to 
the Palatini approach [5(j, it is worth noting that f(R) 
theories were initially investigated in their metric for- 
mulation. In particular, in Ref. \57\. a toy model with 



Since, in this class of models, higher order geometrical terms give 
rise to a quintessential behaviour, it is c usto mary to refer to this 
scenario also as curvature quintessence l23l . 



f(R) = foR n and no matter term was discussed in de- 
tail and confronted with the SNela Hubble diagram and 
the age of the universe. Choosing between the Palatini 
and the metric approach for the f{R) theories is an open 
problem and a definitive answer is likely far to come. 
Actually, a significative advantage of the Palatini for- 
mulation with respect to the metric approach to f(R) 
theories is related to the mathematical simplicity of the 
dynamical equations that turns out to be of second or- 
der. On the other hand, the metric approach leads to a 
fourth order nonlinear differential equation for the scale 
factor a(t) that, in general, cannot be solved analytically 
even for the simplest expressions of f(R). Looking for 
numerical solutions is a difficult task because of the large 
uncertainties on the parameters determining the bound- 
ary conditions. 

It is interesting to wonder whether it is possible to over- 
come these difficulties with a radically different approach 
to f(R) theories in the metric formulation. As we will 
show, the equation governing the dynamic of the universe 
in higher order theories of gravity may be considered as 
a linear third order differential equation for f{R{z)) as 
function of the redshift z. Moreover, the boundary con- 
ditions for f(R(z)) are easily set on the basis of physical 
considerations only and are thus not affected by measure- 
ment errors. Since the Ricci curvature scalar R may be 
also expressed as a function of z, it is then immediate 
to get f(R) and thus determine the gravity Lagrangian. 
The price to pay is to choose an ansatz for H(z), with 
H = a/a the expansion rate. Such an approach makes 
it possible to directly reconstruct f(R) from the data 
since H(z) may be determined, for instance, from the 
SNela Hubble diagram and/or the gas mass fraction in 
galaxy clusters in a model independent way. Rather than 
proposing a theory and testing it against the data, we 
follow here the opposite way around : we start from the 
data and determine which is the f(R) theory that repro- 
duces those data. In this sense, we are facing an inverse 
problem. 

Although it can be conceived as an observationally - 
based method, our approach allows to tell something 
more. As we will see, what is needed to reconstruct f(R) 
is an expression for H(z), whatever is its underlying mo- 
tivation. As a consequence, one could adopt for H(z) 
what is predicted by a given dark energy model (what- 
ever it is) and determine what is the f(R) theory which 
gives rise to the same dynamics (i.e. the same expansion 
rate and scale factor). As a consequence, we will show 
that there is an intrinsic degeneracy among f(R) theo- 
ries and dark energy models so that these two apparently 
radically different explanations of cosmic acceleration are 
reconciled as two different aspects of the same underlying 
physics. 

The outline of the paper is as follows. In Sect. II we 
derive the cosmological equations from a generic f(R) 
theory. Sect. Ill is devoted to a detailed explanation of 
the method to determine f(R) from a given expression 
for the expansion rate H as function of the redshift z. An 
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ansatz for H(z) is suggested starting from astrophysical 
data or from theoretical motivations. The former case 
is discussed in Sect. IV and makes it possible to recon- 
struct f(R) directly from the data. The latter possibility 
instead allows to reconcile f(R) theories with popular 
dark energy models and is discussed in Sect.V. A sum- 
mary of the results and the conclusions are presented in 
Sect. VI. In the Appendix, we comment on the stability 
of f(R) theories in the metric formulation. 



account by introducing an effective fluid and a coupling 
of the matter term to the geometry through the function 
f'(R)- Note that the ansatz f(R) = R + 2 A reproduces 
the standard General Relativity with a cosmological con- 
stant term. 

To get the equations governing the dynamics of the 
universe, we have to choose a metric. The assumed homo- 
geneity and isotropy of the universe motivate the choice 
of the Robertson - Walker metric. The modified Fried- 
mann equations are |U : 



II. CURVATURE QUINTESSENCE 

Although being a cornerstone of modern physics, Ein- 
stein General Relativity has been definitely experimen- 
tally tested only up to the Solar System scales. As such, 
it should not come as a surprise if some modifications 
could be needed on the larger cosmological scales. Fur- 
thermore, any attempt of formulating a quantum theory 
of gravitation claims for a revision of the Einstein theory 
which could only be attained by generalizing the grav- 
ity Lagrangian. Motivated by these considerations, it is 
worth considering a generic fourth order theory of gravity 
given by the action : 



A= / d i xV^g[f(R) + C m ] 



(1) 



with g the trace of the metric, f(R) a generic function of 
the Ricci scalar curvature R and C m the standard matter 
Lagrangian. Hereafter, we adopt units such that 8irG = 
c = 1. The field equations may be obtained by varying 
with respect to the metric components and can be recast 
in the expressive form j^J H3 : 



T 



(to) 



G al3 - R a /3 ~ -Rg a /3 - T^p rv) -r-L a/3 



(2) 



where we have defined a stress - energy tensor for the ef- 
fective curvature fluid : 



T 



(curv) 



a{3 



{9afs [f(R)-Rf'(R)}/2- 



f'(R) 

+ /W" {9*n90» - 9ap9v»)} (3) 

and the matter term enters Eq.J5J through the modified 
stress - energy tensor : 



T, 



(m) 



a/3 



T$/f(R) 



(4) 



with the standard minimally coupled matter stress - 
energy tensor. Here and in the following, we denote with 
a prime the derivative with respect to R and with a dot 
that with respect to cosmic time t. Eqs.J2J - (QJ shows 
that the effect of geometrical terms may be taken into 



H z + — = - 



f'(R) 



2 a - 

a 



H 2 



k 



■Pm) 



(5) 



(6) 



where ait) is the scale factor, H = a/a is the Hubble 
parameter, (p m tPm) are the matter-energy density and 
pressure and we have defined the same quantities for the 
effective curvature fluid as : 



' ( l[f(R)-Rf'(R)}-3HRf"(R)} , 



f'{R) 12 



(7) 



Pcurv — l^curv Pcurv (8) 

with the effective barotropic factor given by : 



= -I 



Rf"{R) + R Rf"'(R) - Hf'(R) 



[f(R) - Rf'(R)} /2 - 3HRf"(R) 



(9) 



Note that the coupling with geometry makes the Fried- 
mann equations nonlinear in the scale factor a(t) since 
this quantity enters the definition of the energy density 
and the pressure of the curvature fluid. Moreover, the 
Hubble parameter H (and thus the scale factor) deter- 
mines how the Ricci scalar evolves with time through the 
following constraint equation : 



f"(R) 



R + 6 



H 2 







R 



6 [H + 2H 2 



(10) 



where we have used the definition of H and assumed 
f"(R) 7^ 0. Note that this latter assumption means that 
we are indeed considering a generalization of Einstein 
theory since f"(R) = implies f(R) oc R that is indeed 
the standard General Relativity. 

There is still another equation that has to be consid- 
ered. Applying the Bianchi identity to Eq.J2J), we get a 
conservation equation for the total energy density : 
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Ptot + 3H(p t ot + Ptot) = . 

Starting from this general scheme, we want to reconstruct 
the form of f(R) from the Hubble parameter as a function 
of the redshift z. According to the present day observa- 
tions 0, 0j @> w e can assume a spatially flat universe 
(k = 0) filled with dust matter. 

Inserting p tot = p curv + p m /f'(R) and modelling the 
matter as dust (i.e. p m — 0), we get : 



3H(1 



)Pc 



1 



J curv JHcurv 



f'(R) 

df'(R) 

Pm dt 



(p m + 3H pr- 



ill) 



Since there is no interaction between the matter and cur- 
vature fluid, we may assume that the matter energy den- 
sity is conserved so that : 



Pm — QmPc 



A I 



(1 + z) 3 



(12) 



It is convenient to change variable from cosmic time t to 
redshift z. To this end, one has simply to use : 



d 
Tt 



a + z)H- 



dz 



With this rule, Ea. H10|) can be written, for k — 0, as : 



R 



-6 



2H 2 



(l + z)H 



dH 

dz 



(15) 



Although straightforward to derive, we find useful to re- 
port the explicit expression of the derivatives of f{R) 
with respect to R as function of the same derivatives 
with respect to z. It is : 



f'(R) 



dR\' x df 



dz 



dz 



UR\ 


~ 2 d 2 f 


fdR\ 


~ 3 d 2 R df 


(dz 


dz 2 


\dz 


dz 2 dz 



(16) 



(17) 



with z = l/a — 1 the redshift (having set a(to) = 1), Qm 
the matter density parameter 2 and hereon quantities la- 
belled with the subscript refers to present day (z = 0) 
values. Inserting Eq. l|12|) into Ea. (|llf> . we get the fol- 
lowing conservation equation for the effective curvature 
fluid: 



f"(R) 



dR\ " 3 d 3 f 
dz J dz 3 

4 

'3 



dR 

dz 



dR 



- 5 f^R 
dz J \ dz 2 
d 2 R d 2 f d 3 R df 



dz 2 dz 2 dz 3 dz 



df_ 

dz 



(18) 



3H(l 



3H 2 n M (i + z y 



Rf"{R) 
[f'{R)f 



(13) 



Actually, the continuity equation and the two Friedmann 
equations are not independent since it is possible to show 
that the dynamics of the universe is completely deter- 
mined by the two cosmological equations [5^. Let us 
then consider only Eqs.© and ©. Using the definition 
of H , we may conveniently reduce the system formed by 
these equations to the following single equation : 



H 



Pi, 



f'{R) 



+ (1 + W curv )p c 



= 



Using Eq.© for 



and Ea. ll2|) for p m , we get : 



H = - 



+R 



{3H 2 n M (l + z) 3 + Rf"(R) + 
Rf'"{R)-Hf"{R)\) . 



2f'(R) 



(14) 



2 Note that Hjvf is defined in terms of the usuai critical energy 
density, but it is not forced to be unity in a spatially flat universe 
since the role of Q\ is now played by Q C urv which can be formally 
derived starting from the curvature quantities. 



Using Eqs.lO- {TBI, it is possible to rewrite Ea. lfHI) 
as a differential equation containing only the unknown 
functions H(z) and f(z) and their derivatives with re- 
spect to z. Note that we are using the abuse of notation 
f(R(z)) = f(z). Some simple but lenghty algebra allows 
to rewrite Ea. p4|) as : 



dz A dz z dz 



with : 



-3H 2 fl M (l + z) 3 
(19) 



Hi = 



dz 



R-RH 



3 ,« 

dz 



dR 

dz 



d 2 R 

dz 2 



2(1 + z)H dH ( dR 
dz V dz 



d 2 R 
~d^ 

-1 



n 2 



R-RH 



dR 

dz 



- 3R 2 



dR 



dR 

dz 



h 3 = r'\-t 

dz 



d 3 R 

dz 3 



- 4 d 2 R 
dz 2 



(20) 



(21) 



(22) 
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The derivatives of R may be evaluated from Eq. (|15|) . For 
instance, differentiating with respect to z, we get : 



dR 

dz 



-(1 + *) 



dH 

dz 



H 



2, dH (l + z ) d2H 
dz dz 2 



and so on for higher order derivatives. For completeness, 
we also report the following useful relations : 



• ,„ . „ dR 
R = -(l + z)H— 

dz 



(24) 



R-RH = 6(1 



H 



z)H 2 |3(l + z) 2 



dH d 2 H 

dz dz 2 
dH] 



dz 



(25) 



Inserting Eqs.© - @ int o Ea. ipi) . we get a cumber- 
some differential equation which we do not report here 
for sake of shortness. Our task is now to solve this equa- 
tion in order to obtain a form of f(R(z)) from the Hubble 
parameter H{z). 



III. SOLVING WITH RESPECT TO f(z) 

In the usual metric approach to f(R) theories in cos- 
mology, one chooses an expression for f(R), uses Ea. l|15|) 
to replace R with H and finally get Ea. i|14|) as a differ- 
ential equation for the scale factor. Unfortunately, here 
the game is over. Indeed, Ea. ll9|) is a fourth order non- 
linear differential equation for a(t) which is difficult to 
solve analytically also for the simplest models of f(R) 
such as power law or logarithmic. To overcome this prob- 
lem, one could impose an analytical ansatz for a(t) and 
look whether and under which conditions the chosen ex- 
pression solves Ea.l|14|). Actually, this approach is quite 
unsatisfactory since there are no hints which may suggest 
a possible expression for a(t) so that one has to perform 
a blind search with the risk of wasting time without end- 
ing with any successful result. Moreover, even if such a 
solution were found, it should be of limited applicabil- 
ity so that drawing any conclusion on the viability of a 
given f(R) theory on the basis of this particular solution 
is dangerous. A possible way to overcome this problem 
is to resort to numerical solutions. However, also this 
approach is plagued by its own problems. Actually, nu- 
merically solving Ea. H14|) demands for the knowledge of 
boundary conditions. Since the equation is of fourth or- 
der, one should give the present day values of the scale 
factor and its derivatives up to the third order. From an 
observational point of view, this means that one should 
set the values of a , the Hubble constant H , the de- 
celeration parameter go and the jerk parameter jo [EjJ 
or one should take into account another set of param- 
eters (see e.g. |6C| for a discussion). While in a flat 



universe we may set oo = 1 and Hq is reasonably well 
constrained, setting qo and jo is quite complicated 3 . It 
is worth stressing that qo and jo should be estimated in 
a model independent way. To this end, one may Tay- 
lor expand the scale factor up to the fourth order (or to 
higher orders) and fit the corresponding luminosity dis- 
tance to the SNela Hubble diagram 62]. While this is 
technically possible, the results are affected by quite large 
errors (especially for jo) and, moreover, are somewhat de- 
pendent on the order to which the series is truncated 4 . 
As a consequence, one should numerically solve Ea. l|14|) 
for a quite large set of boundary conditions thus leading 
to a so great uncertainty on the reconstructed a(t) that 
it remains practically undetermined. 

Given this situation, it is worth asking whether a rad- 
ically different approach is possible. To this end, let us 
first observe that what is indeed determined from the 
data is the Hubble parameter H (z) since this is the quan- 
tity entering the definition of both the luminosity and an- 
gular diameter distances which are tested by the SNela 
Hubble diagram and the data on the gas mass fraction 
in galaxy clusters respectively. Actually, this is what 
is usually done when matching any dark energy model 
to this kind of astrophysical data. Given a background 
cosmological model, one first compute the correspond- 
ing H (z; p) and then constrains the model parameters 
p by comparing with the data. Moreover, H(z) could 
also be directly recovered from the luminosity distance 
in a model independent way |64| although this proce- 
dure is still affected by large errors. Motivated by these 
considerations, we try to go backward from the data to 
f(R) by assuming an ansatz for the Hubble parameter 
H (z) rather than for f(R) itself. Inserting an ansatz for 
H(z) into Eqs.(|2U|)- Eq.jnj) can be seen as a third 
order differential equation for the function f{z). Since 
R(z) may be straightforwardly evaluated from Eq. H15|) 
for a given H(z), we may then obtain f(R) by eliminat- 
ing (numerically) the redshift z from the solution f(z) 
of Ea. (|19|l . Note that considering f(z) as unknown has 
two immediate advantages. First, the equation is of third 
rather than fourth order so that the numerical solution is 
easier to find. Second, the equation is linear in f(z) and 
its derivatives so that, for given boundary conditions, the 
solution exists and is unique. 

To numerically solve the differential equation for f(z), 
one has to set the boundary conditions, i.e. the values 
of / and its first and second derivatives with respect to 
z evaluated at z = 0. Here lies another advantage of this 
approach : the boundary conditions may be chosen on 
the basis of physical considerations only. To this end, let 
us first remind that it has been shown that, in order to 



3 Note that the boundary conditions may also be set using the 
statefinder parameter ro 10 ll instead of the jerk jo- 

4 To be more precise, the best fit values of go an d jo are almost 
unaltered, but the uncertainties get larger as the order of the 
series expansion increases. 



7 



not contradict Solar System tests of gravity, a whatever 
f(R) theory must fulfill the condition f"{R ) = [g|. 
Using Eq. fTTjl , we thus get the constraint : 



f"(Ro) 



(dR\ 


~ 2 d 2 f 


UR\ 


3 d 2 R df 


[dz 


dz 2 


[dz 


dz 2 dz 



J z=0 



A second condition may be obtained considering Eq.JjJJ. 
Let us rewrite it introducing back the coupling factor 
8ttG and k = 0. It is : 



H 2 



8ttG 



/'OR) 



This equation shows that the coupling of the matter with 
geometry through the function f'(R) is equivalent to re- 
define the Newton gravitational constant G as G/f'(R) 
so that we get the well known result that in f(R) the- 
ories G is redshift (and hence time) dependent. Indeed, 
at z = 0, the effective gravitational constant G/f'(Ro) 
must be equal to G so that we get : 



f'(Ro) 



dR 

dz 



- l d[_ 

dz 



= 1 



J z=0 



Combinining the constraints on f'(Ro) and f"(Ro), we 
get the following boundary conditions : 



z=0 



dR 

dz 



(26) 



z=0 



Evaluating Eq.JSJ) at z 
f(Ro), we finally get: 



and solving with respect to 



f(z = Q) = f(Ro) = 6H 2 (l - M ) + Ro 



(28) 



Summarizing, to reconstruct f(R) from the data, we 
adopt the procedure schematically sketched below. 

1. Assume an expression for H{z) and determine the 
model parameters from fitting to the data. 

2. Compute R(z) from Eq. (fL5|) . 

3. Use Eqs.JU- (E3 to write Eq.© for f(z). 

4. Solve numerically the above equation with the 
boundary conditions given by Eas. i|26|) - l|28|l . 

5. With the aid of the derived expressions for R(z) and 
f{z), eliminate z from f(z) and finally get f(R). 

This quick pipeline 5 makes it possible to obtain an ex- 
pression for f(R) that is observationally well founded 
since the f(R) theory thus reconstructed fits the same 
dataset used to determine the parameters entering H[z). 
It is worth stressing that such a result has been obtained 
without the need to solve Ea. (|14|) with respect to the 
scale factor a(t). Moreover, any arbitraryness in the 
choice of an expression for f{R) has been removed. 



IV. DETERMINING f(R) FROM THE DATA 



dz 2 



d 2 R 
Ik 2 



(27) 



A comment is in order here. We have derived the present 
day values of df j dz and d 2 f /dz 2 by imposing the consis- 
tency of the reconstructed f(R) theory with local Solar 
System test. One could wonder whether tests on local 
scales could be used to set the boundary conditions for a 
cosmological problem. It is easy to see that this is indeed 
meaningful. Actually, the isotropy and homogeneity of 
the universe ensure that the present day value of a what- 
ever cosmological quantity does not depend on where the 
observer is. As a consequence, an hypothetical observer 
living in the Andromeda galaxy and testing gravity in his 
planetary system should get the same results. As such, 
the present day values of df /dz and d 2 f /dz 2 adopted by 
this hypothetical observer are the same as those we have 
used based on our Solar System experiments. Therefore, 
there is no systematic error induced by our method of 
setting the boundary conditions. 

Finally, to set f(z = 0), let us evaluate the present day 
value of Pcurv using the definition Q and Eas. l|26|l and 
i P7| l. We get: 



,(z = 0) 



f(Ro) - Ro 



A key ingredient in the procedure outlined above is 
the choice of an analytical expression for the dependence 
of the Hubble parameter on the redshift. Rather than 
choosing a somewhat motivated ansatz, one should resort 
to an empirical determination of H{z) from the data. 
When used as input for our pipeline, this function makes 
it possible to reconstruct f(R) directly from the data thus 
avoiding any systematic bias or theoretical prejudice. We 
will consider here two different approaches to implement 
an observations - based determination of f(R). 



A. f(R) from D L (z) 

Fitting to the SNela Hubble diagram is nowaday a 
standard tool to investigate the viability of a given cos- 
mological model. The essence of the method relies on the 
well known relation : 



H(z) = 51og£> L (z) + 25 



(29) 



5 A fast and efficient code implementing this procedure, written 
for the Mathematica 4-1 software, is available on request. 
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being [i the distance modulus of an object at redshift z 
and Dl its luminosity distance (in Mpc) which is deter- 
mined by the Hubble parameter as follows 6 : 



D L (z) = (l + z) 



la H(C) • 

This equation may be inverted to give |63j | : 



(30) 



A R = 



382 <Jo 



(l + z)5i-So 



3(1 + z)S 2 ai 



5 - (1 + z)5 l 



Inverting Ea. (|29J) . we may express (do, 61,62) as functions 
of the distance modulus fi(z) and its derivatives up to the 
second order. Following the same procedure for propa- 
gating the errors, we get : 



H(z) 



d 

dz 



Dl(z) 



1 



(31) 



so that H(z) may be determined from Dl(z). Because 
of Eq. l|29fl . the luminosity distance may be empirically 
determined by measuring the distance modulus for a class 
of standard candles thus allowing a direct reconstruction 
of H (z) from the data in a model independent way. 

It is worth wondering how much accurate the determi- 
nation of (j,(z) must be in order to efficiently recover H(z) 
and thus f(R). A quantitative answer requires detailed 
simulations and it is outside the scope of the present pa- 
per. Actually, we may get a qualitative understanding of 
the problem by considering the error on R(z) due to the 
measurement uncertainties on n(z). To this end, let us 
first insert Eq. into Eq. (JTSJl to get : 



R(z) = 



6(1 



[D L (z)-(l + z)D' L (z)Y 



(32) 



where, only in this section, the prime denotes the deriva- 
tive with respect to z. If the errors on Dl, D' l and D'[ are 
Gaussian distributed, we can determine the uncertainty 
<tr on R by the usual rules. This is likely to be not the 
case. Nonetheless, propagating the errors on Dl, D' l and 
D'l should give us an order of magnitude to estimate an 
that is enough for our aims. To this end, let us denote 
with (<5o ? ^1 j <^2 ) the quantities (Dl, D' L , D'[) respectively 
and let (ctq, 01, 02) the corresponding uncertainties. A 
naive estimate of ctr may be obtained as : 



or = 



dR 


2 


dR 


2 


dR 


d8o 


a 2 + 


d5x 




d6 2 



Inserting Eq. lIML'l) into the above relation, we get : 



or = 



6(l + z) e 



[6 Q -(l + z)8 1 



X A; 



(33) 



with : 



00 = k 6 cr p 



(7i = kJii{a l o + <5o°"mi 



(35) 



(36) 



Remember that we are using units such that c = 1. 



a 2 = k^(k^\ + M2 )V 2 + k(kaf n + a 2 2 ) 2 S 2 , (37) 

having set k = In (10)/5, /ii = d/jt/dz, /12 = d 2 [ijdz 1 and 
denoted with (cr„ , , a^ 2 ) the measurement uncertain- 
ties on (/i, (ii, 112) respectively. Note that Eqs.JSSJl - IJ57JI 
clearly show that the errors on the luminosity distance 
and its derivatives are correlated, while Eq. i|33[l has been 
obtained considering the errors as uncorrelated. As a 
consequence, Eq. l|33|l underestimates <7r. Nonetheless, 
some interesting conclusions on an may be drawn. First, 
we note that or oc (1 + z) 6 so that the error on the recon- 
structed Ricci scalar quickly increases with the redshift. 
This conclusion is further strengthned when considering 
that (80,81,82) become larger and larger with Dl- Al- 
though a detailed investigation (with the aid of simula- 
tions) is needed, one should argue that the distance mod- 
ulus must be measured with a tiny uncertainty in order 
to reduce as more as possible <tr. Moreover, Eas. H32|) and 
(|33|l show that both R and <jr depend on the derivatives 
up to the second order of the distance modulus with re- 
spect to the redshift. Actually, observationally constrain- 
ing d\ijdz and d 2 H/dz 2 is a quite complicated task likely 
demanding for very large samples of SNela. Since recon- 
structing f(R) needs first a determination of R(z), we 
may thus conclude that a completely model independent 
determination of f(R) from the luminosity distance data 
is plagued by too large uncertainties and is far to come. 

A possible way to escape these problems is to resort 
to some parametrized expression of Dl(z) so that one 
has only to constrain a limited set of quantities rather 
than recover a function (i.e., determine an infinite num- 
ber of unknowns). To this end, a series expansion |64j 
or a more versatile fitting function [6(j have been pro- 
posed. Inserting the expression for Dl into Eas. l|32|l and 
(|33|l . one could evaluate both R and gr as function of 
the parameters assigning the luminosity distance so that 
the problem of reconstructing R(z) reduces to the de- 
termination of these parameters. Similar considerations 
also hold for ,f(z). Calibrated simulations are needed to 
investigate the viability of this promising approach. 
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B. Polynomial fit to the dark energy density 

The problem to recover f(R) from the luminosity dis- 
tance in a model independent way is partially related to 
the need of determining from the data not only Dl(z), 
but its derivatives too. To overcome this difficulty, one 
should resort to a less ambitious program choosing a 
parametric ansatz for H(z) which should be as general 
as possible. A useful expression is : 



Ax + A 2 x + A 3 x 2 



(38) 



with x = 1 + z and (Ai, A%, A 3 ) parameters to be deter- 
mined from the data. Note that the case [A\,A2,A%) = 
(0,0,0) corresponds to a flat matter only universe, while 
the case A x ^ and (A 2 ,A 3 ) = (0,0) gives H(z) for the 
ACDM concordance model. In general, Ea. (|38[) refers to 
a universe containing a matter term and a second compo- 
nent whose energy density is approximated by a second 
order polynomial fit. 

This expression has the remarkable property of sim- 
plicity. Actually, by inserting Ea. i|38|) into Ea. l|32() . we 
get: 



R = -3iJ 2 (4A! + 3A 2 + 2A 3 + Sl u ) , (39) 
while the boundary conditions turn out to be given as : 

f(z = 0) = -3H 2 (4A X + 3A 2 + 2A 3 + 30 M - 2) , (40) 



df/dz\ z=Q = -3H*(3A 2 + 4A 3 + 3V M ) , 



(41) 



d 2 f/dz% =0 = -3H 2 (2A 3 + 3fl M ) . (42) 

Denoting with <jj the uncertainty on A4 and with 
(a Hi cm) those on Hq and respectively, a naive esti- 
mate of the error on R is given by : 



a R = 3H 2 ^ (12/ H 2 )a 2 H + 16a 2 + 9a 2 + 4<r| + a 2 M 

having assumed that the errors are statistical and un- 
corrected. It is worth noting that <tr turn out to 
be independent of the redshift z so that a reliable re- 
construction of the evolution of the scalar curvature 
only demands for lowering the errors on the parameters 
(Ho, Om) Ai, A 2 , A 3 ). This could be achieved by increas- 
ing the sample of SNela or reducing the measurement 
and systematic errors on the distance modulus or by a 
judicious use of priors on (Hq, Q.m)- 

Having reconstructed R(z), one has only to determine 
f(z). Inserting Eq.{3HJ| into Eq.JTJJJ, we get : 



+ V 2 (z)^ + Vi(z)^ = -3H 2 n M (l + zf 

. . (44) 

where Vi are cumbersome functions of the redshift z and 
the parameters (Qm, A±, A%, A3) that may be computed 
inserting Eq.^SJ mto Eas.p0" )l - l|23 |l . We do not report 
them for sake of shortness. Ea.(|44|) may be straightfor- 
wardly solved numerically using the boundary conditions 
(f4*U|) - 1|4^|) so that f(z) may be obtained and coupled with 
the above reconstructed R(z) to finally determine f(R) 
from the data. Estimating what is the error on f(R) is 
quite complicated and it is likely that a numerical anal- 
ysis is needed. 



V. f(R) THEORIES AND DARK ENERGY 
MODELS 

As discussed in the previous section, H(z) may be 
directly reconstructed from the data on the luminosity 
distance. Although in principle possible, this method 
is nowaday still not feasible since it is affected by quite 
large errors. Furthermore, H(z) is thus recovered only 
on the redshift range probed by the data used so that 
a potentially dangerous extrapolation is needed to go to 
higher z. This problem forces us to adopt for H(z) a the- 
oretically rather than observationally motivated ansatz. 
To this end, it is worth referring to the literature where 
the accelerated expansion of the universe is usually ex- 
plained proposing a wide variety of dark energy mod- 
els. Although different in their physical background, all 
these models share the property of well fitting the same 
dataset so that, from this point of view, we could adopt 
for H(z) the expression corresponding to one of these 
models. This choice also allows us to stress one interest- 
ing point. Since the procedure above sketched makes it 
possible to recover f(R) from H(z), we are able to con- 
struct an f(R) theory of gravity which gives the same 
dynamics (i.e., the same expansion rate and scale fac- 
tor) of a whatever dark energy model. We will show this 
explicitly for two popular dark energy models, namely 
quintessence with constant equation of state referred to 
as quiessence and the Chaplygin gas. 

As a preliminary remark, it is worth discussing what 
are the dimensions of the different quantities involved 
in the units (with 8nG = c = 1) we are adopting. To 
this end, let us first consider Ea. (|10fl . Since H is ex- 
pressed in s _1 , the Ricci scalar turns out to be mea- 
sured in s~ 2 . From Eq.© and the consideration that 
Pcurv and p rn have the same dimensions, we conclude 
that f'(R) — df/dR is dimensionless and hence / has 
the same dimension as R. Finally, the energy density is 
expressed in s~~ 2 as can be inferred from the expression 
of the critical density that is p cr it = 3H 2 in these units. 
It is straightforward to check that this is consistent with 
what is obtained for p curv from Eq.lQ). Hereon, we mea- 
sure time in units of 1/Hq so that Hq = 1 and all the 
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FIG. 1: Reconstructed f(z) for quiessence models with Qm = 
0.3 and three different values of the dark energy barotropic 
factor, namely w = —0.5 (short dashed), w = — 1 (solid) and 
w = —1.5 (long dashed). We report If = In ( — /) rather than 
/ to better show the results. 
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FIG. 2: Reconstructed f(R) for quiessence models with Qm = 
0.3 and three different values of the dark energy barotropic 
factor, namely w — —0.5 (short dashed), w — —1 (solid) and 
w = —1.5 (long dashed). We report If — ln( — /) as function 
of IR = In (-R) rather than f(R) to better show the results. 



quantities we are interested in, namely R and f(R), are 
dimensionless that is a useful feature when dealing with 
numerically solving differential equations. 



A. Quiessence from f(R) 

As a first straightforward application, let us consider 
the ansatz : 



H{z) = H ^Jn M (l + zf + n x (l + z)3d+-) (45) 

with fix = (1 — Hjvr) and w a constant parameter. 
Ea. (|45|l gives the Hubble parameter for the so called 
quiessence models (or QCDM) where the acceleration of 
the universe is due to a negative pressure fluid with con- 
stant barotropic factor w. This is the easiest generaliza- 
tion of the cosmological constant which is obtained for 
w = —1. Quiessence has been successfully tested against 
the SNela Hubble diagram and the CMBR anisotropy 
spectrum that has made it possible to severely constraint 
the barotropic factor w (6^ ■ It is interesting to note that 
these constraints extend into the region w < — 1 so that 
models violating the weak energy condition are allowed. 
This class of models, dubbed phantom models, are af- 
fected by serious problems with the growth of perturba- 
tions and are therefore worrisome. 

Inserting Eq. (|45|l into Eq. i|15|) , we get : 



R = -3ff 2 \q m (1 + zf +Sl x (l- 3w)(l + zf il+w ^ 

(46) 

Note that R is always negative as a consequence of the 
signature {+,—,—,—} adopted. If we had used the op- 
posite signature, Ea, 14t)|) is the same, but with an overall 
positive sign. The ansatz in Eq. I|45|l leads to the following 
equation for f{z) : 



d 3 f 



d 2 f 



df 



dz 



dz 



dz 



where Qi{z) may be evaluated by inserting Ea. l|45|) into 
Eas. l|20|l - i|25|) . The resulting cumbersome functions of 
the redshift z and the model parameters (Om , w) are 
not reported here for sake of shortness. Eq. (|47|l may be 
easily solved numerically using the boundary conditions 
|[11)J-JSBJ|. The result is shown in Fig.[T]for models with 
fW = 0.3 and three different choices of the barotropic 
factor w. Note that we have plotted —f(z) rather than 
f(z) since / turns out to be negative because of the sig- 
nature adopted. Moreover, we use a logarithmic scale to 
better show the results. For z < 0.5 the three curves 
show a different behaviour, while, for higher z, the shape 
of f(z) is unaffected by w which only acts as a scaling pa- 
rameter. This is also true for Qm ■ changing this param- 
eter only shifts up or down the curves plotted in Fig.^ 
Inverting numerically Ea. (|46|l . we may obtain z — z(R) 
and finally get f(R) shown in Fig.|2for the same models 
considered above. It turns out that f(R) is the same for 
different models for low values of R and hence of z. This 
is a consequence of the well known degeneracy among 
different quiessence models at low z that, in the stan- 
dard analysis, leads to large uncertainties on w. Because 
of this degeneracy, models with different barotropic fac- 
tors may not be discriminated and are therefore dynam- 
ically equivalent. This is reflected in the shape of the 
reconstructed f(R) that is almost w - independent in this 
redshift range. 

Fig.|21may suggest that ln(— /) is a quasi- linear func- 
tion of In (— R) that is not the case. Actually, we have 
checked that the following empirical function 



h [In (-R)] h [1 + In (-R)] h + h (48) 
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approximates very well the numerical solution provided 
that the parameters I2, 13, 14) are suitably chosen for 
a given value of w. For instance, forio = -l (the cosmo- 
logical constant 7 ) it is : 

[hMMM) = (2.6693,0.5950,0.0719,-3.0099) . 

The error in this case is of the order ~ 2 — 10% for z < 1, 
but remains smaller than ~ 4% up to z — 10. We have 
checked that Eq. (|48|) works very well also for other values 
of w in the range (—1.6,-0.6) with an error 8 depending 
on z and w, but always smaller than ~ 10%. 

Some comments are in order as final remarks. First, 
remember that we are using units with 8ttG = c = 1 
and 1/Hq = 1. As a consequence, the values of the fit- 
ting parameters should be adjusted if physical units are 
used. However, this does not change the shape of the ap- 
proximating function since the numerical solution does 
not depend on the adopted units. Let us also stress that 
Eci. (|48|) is only a fitting function tested over the large, 
but still limited redshift range (0,10). Extrapolation to 
higher z may introduce large errors and should therefore 
be avoided. 



B. Chaplygin gas as f(R) - model 

Recently, much attention has been devoted to models 
where a single fluid with an exotic equation of state ac- 
counts for both dark matter and dark energy. Usually 
referred to as unified dark energy models, such models 
have the nice feature to solve two problems in a single 
step and have thus attracted a lot of interest. As a pro- 
totype example, we consider the generalized Chaplygin 
gas (GCG) whose equation of state is : 



P = -A/ P ° 



(49) 



with A and a positive parameters to be determined from 
the data. Assuming a spatially flat universe, the Hubble 
parameter is : 



H(z) = H [As + (1 - A s )(l + z) 3il+ ^\ 2(1+Q) (50) 

with A s = A/p(z = 0). Note that A s — 1 reduces to a 
universe with the cosmological constant as unique com- 
ponent, while A s = describes a matter dominated uni- 
verse. Motivated by this consideration, we will only take 



7 The case of the cosmological constant may be solved analyti- 
cally giving f(R) = -R + 2A which is quite different from Eg. 1481 . 
Actually, Eg. 1481 is only an approximating function chosen be- 
cause of its versatilty. Indeed, for vj = — 1, the fitting parameters 
renders the approximating function as similar as possible to the 
exact expression. 

8 Tables with the values of the parameters and the approximation 
error for different w and z are available on request. 



into account models with A s in the range (0, 1) although, 
in principle, nothing prevents A s to be larger than 1. 

We may apply our procedure to recover the f(R) the- 
ory which reproduces the same dynamics of the GCG. 
For simplicity, we concentrate our attention only to the 
case a — 1 that is the originally proposed Chaplygin gas. 
The Ricci scalar turns out to be : 



R 



C(z)- 



3(1- A s )(l + z) 



61 



C(z) 



with : 



C(z) = ^A s + {l-A s )(l + zf 



(51) 



(52) 



Proceeding as for the quiessence case, we get the follow- 
ing equation for f(z) : 



'dz 3 



C 3 (^)^+C 2 (z)^+C 1 (z) 



] d? 



dz 



(53) 



where Cj(z) are functions of the redshift z and A s 
which are derived by inserting Ea. l|50|l into Eas. (|20|) - 
(|25|l . The numerical solution of Ea. (|53|) is plotted in 
Fig-El while Fig.0] shows the corresponding f(R) ob- 
tained by eliminating z using Eq. l|51|l . An important 
caveat is in order here. Although ft^j is identically 1 in 
the usual approach to Chaplygin gas, one has to choose 
a value of 7^ 1 since the Friedmann equations have 
been modified. Moreover, this choice should be based 
on model - independent determination of this parameter 
which also enters the determination of the boundary con- 
dition f(z = 0) through Eq.®. We set Q, M = 0.3 as for 
the best fit concordance ACDM model since this value is 
also in agreement with the model independent estimates 
coming from the abundance of galaxy clusters. 

The similarity among the curves for different values 
of A s is striking, but not fully unexpected. Actually, at 
larger z (where the curves are almost overlapping in these 
logarithmic plots), the Chaplygin gas reduces to a matter 
only universe whatever is the value of A s so that the 
dependence on this parameter is washed out. A stronger 
dependence on a may be expected in the GCG scenario 
since this parameter controls the rate of transition from 
a A dominated to a matter dominated universe. 

Finally, we have checked that the numerical solution 
for f(R) is very well approximated by Eq.(@2J with : 

(h,h,h,k) = (1.9814,0.5558,0.2665,-2.5337) 

for the model with (Q M ,A S ) = (0.3,0.75). We stress 
again that Ea. l|48f) is only an empirical fitting function 
so that it is dangerous to draw any physical implica- 
tions from the fact that the same functional expression 
for f(R) works well for both quiessence and Chaplygin 
gas models. Actually, such a result could be somewhat 
explained noting that both classes of models are designed 
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FIG. 3: Reconstructed f(z) for Chaplygin gas with Om = 
0.3 and three different values of the normalization constant, 
namely A s = 0.25 (short dashed), A a — 0.50 (solid) and A 3 — 
0.75 (long dashed). We report If — ln( — /) rather than / to 
better show the results. 

to fit the same dataset over the same redshift range. 
As such, there is a sort of degeneracy among quiessence 
and Chaplygin gas that could be the reason why Eq. 
works well for both class of dark energy models. 

VI. DISCUSSION AND CONCLUSIONS 

The observed cosmic acceleration may be seen as the 
first signal of a breakdown of the Einstein theory of Gen- 
eral Relativity. Indeed, several theoretical motivations 
claims for modifications of General Relativity. Moti- 
vated by these considerations, much attention has been 
recently devoted to higher order theories of gravitation 
which are obtained by replacing the Ricci scalar R with a 
generic function f{R) in the gravity Lagrangian. Usually 
referred to as f(R) theories, these scenarios make it pos- 
sible to explain the observed cosmic acceleration without 
introducing any scalar field or changing the properties of 
matter term (curvature quintessence). Although phys- 
ically well motivated and mathematically elegant, this 
approach has its own problems. The metric formulation 
of f{R) theories leads to a fourth order nonlinear dif- 
ferential equation for the scale factor a(t) which cannot 
be analytically solved in general even for some simple 
choices of the function f(R). Moreover, a numerical so- 
lution, although possible, is difficult to handle because of 
the large uncertainties on the parameters (do, Hq, go, Jo) 
which determine the boundary conditions. 

It is worth noting that the same equation may also be 
seen as a linear third order differential equation for the 
function f(z). For a given H(z), this equation is easier to 
solve numerically since the boundary conditions may be 
set on the basis of physical considerations only. We have 
therefore developed a straightforward procedure which 
makes it possible to solve for f(R), given an expression of 
the Hubble parameter H(z) and a physically motivated 
choice of the boundary conditions. We stress that our 
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FIG. 4: Reconstructed f(R) for Chaplygin gas with CIm = 
0.3 and three different values of the normalization constant, 
namely A s = 0.25 (short dashed), A 3 = 0.50 (solid) and A s = 
0.75 (long dashed). We report If = ln( — /) as function of 
IR = In (— R) rather than f(R) to better show the results. 



approach to f(R) theories is reversed with respect to the 
usual one. Rather than arbitrary choosing the function 
f(R) and then compare the corresponding model to the 
observational data, we reconstruct f(R) from a quantity 
which (at least, in principle) may be directly estimated 
from the data. As a consequent result, the f(R) theory 
so obtained will intrinsically fit the data and thus it is a 
viable candidate to solve the dark energy puzzle. 

Our method to determine f(R) needs a precise knowl- 
edge of the Hubble parameter H (z). This quantity could 
be derived, for instance, from the luminosity distance of 
SNela, but the uncertainties are still too large thus pre- 
venting a full reconstruction of f(R) directly from the 
data. An easy way to overcome this problem is resorting 
to a parametrized expression of H(z) to be fitted against 
the available data in order to determine its parameters. 
This function may then be used as input for the pipeline 
we have devised. We have explored one possible choice 
also evaluating (the order of magnitude of) the uncer- 
tainty on R(z). Exploring in detail this approach, how- 
ever, demands for simulated datasets in order to estimate 
the error on the reconstructed f(R) and the properties 
of the sample needed for a reliable determination of the 
gravity Lagrangian directly from the data. This further 
issue will be investigated in a forthcoming paper. 

Although implemented as an observations - based 
method, it is possible to use our procedure as a simple 
way to investigate intrinsic degeneracies among popular 
dark energy models and f(R) theories. We have indeed 
shown that from the point of view of dynamics of the 
late time universe, higher order theories of gravity may 
mimic several dark energy models. Actually, using as in- 
put for our pipeline the Hubble parameter corresponding 
to a given model, we may found the f(R) theory which 
reproduces the same dynamics and thus fit the data with 
the same goodness as the given dark energy model. We 
have explicitly shown this for quiessence (dark energy 
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with constant equation of state) and Chaplygin gas (as 
a prototype of UDE models). It is worth noting that 
the same functional expression fit well the numerically 
reconstructed f(R) for both class of models. Although 
unexpected, this result could also be a consequence of the 
fact that both models fit the same dataset so that they 
are forced to resemble each other over the redshift range 
probed by the data. However, the corresponding f(R) 
theories extend this similarity over the full evolutionary 
history of the universe since it is likely that the functional 
shape of f(R) does not change with z even if R evolves. 
The successful application of the method to two radically 
different dark energy models is perhaps the most impor- 
tant result of this paper since it suggests that the crowded 
zoo of dark energy models may be seen as the result of 
our ignorance of what is the correct expression for f{R) 
that has to be inserted in the gravity Lagrangian. Ac- 
tually, observations tell us that something is wrong with 
the old standard picture of the universe, but cannot in- 
dicate what it is. It is then somewhat a matter of taste 
to choose whether we are lacking an ingredient (such as 
a scalar field), or the matter equation of state need to be 
modified (as in UDE models) or rather it is the underly- 
ing theory of gravity which has to be generalized. What 
we have shown here is that these three philosophycally 
different approaches may indeed be considered as distinct 
manifestations of the same physics of the late universe. 

Discovering what is this physics remains an open prob- 
lem demanding for both theoretical investigations and 
observational evidences. From an observer's point of 
view, CMBR anisotropy and polarization spectra and 
the clustering properties of the large scale distribution 
of galaxies are ideal tools to discriminate among the ri- 
val possibilities outlined above. Indeed, while the SNela 
Hubble diagram and the data on the gas mass fraction in 
galaxy clusters probe only the Hubble parameter, both 
CMBR and large scale structure depend on how the evo- 
lution of perturbations takes place in the background cos- 
mological model. Actually, the perturbation equations 
for scalar field quintessence and UDE models are rad- 
ically different from those in f(R) theories even if the 
dynamics of the model (controlled by the Hubble param- 
eter) is the same. From a theoretician's point of view, 
choosing among higher order gravity theories and dark 
energy models is quite difficult because of the intrinsic 
degeneracy we have shown. Nonetheless, constraints on 
f(R) theories may be obtained by considering the full 
evolutionary history of the universe. Since the functional 
expression of the gravity Lagrangian does not change 
with z, a given f(R) theory should not only describe the 
late time universe, but also give rise to an inflationary pe- 
riod. If we will be able to reconstruct f(R) directly from 
the data without resorting to a parametrized expression, 
we could study the picture of the universe assigned by 
this model and discard such f(R) theories which are not 
able to reproduce inflation. 

The procedure we have presented relies on the field 
equations @ that have been obtained varying the La- 



grangian with respect to the metric only. Actually, as we 
have discussed above, f(R) theories may also be stud- 
ied using the Palatini approach where the variation is 
performed with respect to the metric and the connec- 
tion considered as independent variables. Being the field 
equations different, the procedure we have implemented 
does not hold in the Palatini formulation of f(R) the- 
ories. We note, however, that it is still not clear what 
is the correct approach to higher order gravity theories. 
Indeed, the possibility of mimicking dark energy models 
offered by the metric formulation of f(R) theories could 
be considered as an encouraging evidence favouring this 
strategy because of the capability to reproduce all the 
successful dark energy phenomenology. A similar method 
has to be developed also for the Palatini formulation, as 
we are going to do in [69j. 

We would like to conclude with a general comment. 
The method we have presented makes it possible to rec- 
oncile dark energy models and f(R) theories as two dif- 
ferent faces of the same medal. Both rival theories are 
able to reproduce the same dataset and thus to describe 
the same dynamics of the late time universe since we have 
shown that they are distinct manifestations of the same 
underlying physics. Discriminating among them is only 
possible going back into the past to the period of struc- 
ture formation and still before, up to inflation. In our 
opinion, a unified study of both the early and late uni- 
verse is the unique way to understand what is the correct 
physics. Higher order gravity theories naturally offer this 
possibility and may thus be used as powerful light to look 
into the dark sector of the universe. 



APPENDIX: STABILITY OF f(R) THEORIES 

Soon after the first proposals of f(R) theories as al- 
ternative explanations of the cosmic acceleration, several 
criticisms were raised on the stability of such models in 
the metric formulation. In particular, in Ref. [70j. it was 
pointed out that the model f(R) = R + fJ. A /R suffers 
violent instabilities and was argued that something sim- 
ilar should hold in every theory of gravity which leads 
to higher order differential equations for the scale fac- 
tor. This argument has been often used as an evidence 
against the metric formulation of f(R) theories thus mo- 
tivating the great interest dedicated to the Palatini ap- 
proach which avoids this problem. Actually things are 
different. It is, indeed, possible to show that a leading 
role in the determination of the stability of the theory is 
played by the following potential (see, e.g., |jl|) : 



U(Ro) 



f^(Ro) f^ 3 \Ro) 



V a R V a R Q 



f(V(Ro) /( 2 )(i? ) 2 
Ro 2fW(R )fW(R )R 
' 3 3/< 2 >(ifc) 2 
/«(i? ) , 2f(R a )f^(R ) 



3/( 2 >(i?o) 



3/( 2 )(i?o) : 
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Rof (3) (Ro) 
3/( 2 ) (i? ) 2 



(A.l) 



where /W = d % f jdR % and Rq is the solution of the un- 
perturbed field equations. Note that, to be coherent with 
Ref. [x3> here we are using the signature { — , +, +, +} so 
that Ro is positive and the trace of the matter stress - 
energy tensor is negative. Ea. (|A.lJI has been obtained by 
setting R = R + R\ with |i?i| << |i? | and developing 
the equation at the first perturbative order. If U{Rq) 
is positive, since DR\ ~ —d^Ri, the perturbation Ri 
becomes exponentially large and the system is unstable. 
Assuming that the matter is uniformly distributed, we 
may simplify Ea. (|A.ll) setting V Q i?o = and then study 
the sign of U(Ro) for a given f(R) theory. 

For the model f(R) = R+[i 4 /R, U(Rq) turns out to be 
positive so that the theory is indeed unstable. But this is 
not a ge neral result. For instance, the choice f(R) = PR" 
0, EE E3 g^es: 



U(R ) 



(n-2)(2/3Eg- 1 -l) 
3/3n(n - l)i$ 



(A.2) 



Assuming that the coupling constant /? is positive, the 
constraint U(Rq) < reduces to: 



2/3R™- 1 - 1 > for n < and 1 < n < 2 



2/3i?™~ 1 - 1 < for < n < 1 and n > 2 



so that the stability of the theory depends on n and (3. A 
similar analysis can be conducted for the approximating 
function (|48|1 in order to determine the stability of the 
reconstructed f(R). It is thus worth stressing, as final 
remark, that the usual criticism against the metric ap- 
proach to f(R) theories about the stability arguments has 
to be reconsidered on a case by case basis. Therefore, we 
are confident that the procedure we have implemented to 
reconstruct f(R) is meaningful and not affected by any 
systematic problem related to the metric formulation of 
f(R) gravity. 



[1] A. G. Riess et al., AJ, 116, 1009, 1998; S. Perlmutter et 
al., ApJ, 517, 565, 1999; R.A. Knop et al., ApJ, 598, 102, 
2003; J.L. Tonry et al., ApJ, 594, 1, 2003; B.J. Barris et 
al., ApJ, 602, 571, 2004; A.G. Riess et al., ApJ, 607, 665, 
2004 

[2] P. de Bernardis et al., Nature, 404, 955, 2000; R. Stompor 
et al., ApJ, 561, L7, 2001; C.B. Netterfield et al., ApJ, 
571, 604, 2002 D.N. Spergel et al. ApJS, 148, 175, 2003 
R. Rebolo et al., MNRAS, 353, 747, 2004 

[3] S. Dodelson et al., ApJ, 572, 140, 2002; W.J. Percival et 
al., MNRAS, 337, 1068, 2002; A.S. Szalay et al., ApJ, 
591, 1, 2003; E. Hawkins et al., MNRAS, 346, 78, 2003; 
A.C. Pope et al., ApJ, 607, 655, 2004 

[4] R.A.C. Croft et al., ApJ, 495, 44, 1998; P. McDonald et 
al., astro -ph/0405013, 2004 

[5] S. Sasaki, PASJ, 48, L119, 1996; U. Pen, New Ast., 
2, 309, 1997; S.W. Allen, R.W. Schmidt, A.C. Fabian, 
MNRAS, 334, Lll, 2002; S.W. Allen, R.W. Schmidt, S. 
Bridle, MNRAS, 346, 593, 2003; M. Sereno, A&A, 412, 
341, 2003; S.W. Allen, R.W. Schmidt, H. Ebeling, A.C. 
Fabian, L. van Speybrock, MNRAS, 353, 457, 2004 

[6] S.M. Carroll, W.H. Press, E.L. Turner, ARAA, 30, 499, 
1992; V. Sahni, A. Starobinski, Int. J. Mod. Phys. D, 9, 
373, 2000 

[7] M. Tegmark et al., Phys. Rev. D, 69, 103501, 2003; U. 

Seljak et al., astro -ph/0407372, 2004 
[8] B. Rathra, P.J.E. Peebles, Phys. Rev. D, 37, 3406, 1988; 

C. Wetterich, Nucl. Phys. B, 302, 645, 1988 
[9] R.R. Caldwell, R. Dave, and P.J. Steinhardt, Phys. Rev. 
Lett. 80, 1582, 1998. 
[10] P.J.E. Peebles, B. Rathra, Rev. Mod. Phys., 75, 559, 
2003; T. Padmanabhan, Phys. Rept., 380, 235, 2003; V. 
Sahni, astro - ph/0403324, 2004 
[11] JW. Moffat, astro -ph/0403266, 2004 
[12] S. Capozziello, V.F. Cardone, S. Carloni, A. Troisi, Phys. 



Lett. A 326, 292, 2004 

[13] A. Kamenshchik, U. Moschella, V. Pasquier, Phys. Lett. 
B, 511, 265, 2001; N. Bilic, C.B. Tupper, R.D. Viollier, 
Phys. Lett. B, 535, 17, 2002; M.C. Bento, O. Bertolami, 
A. A. Sen, Phys. Rev. D, 67, 063003 

[14] G.W. Gibbons, Phys. Lett. B, 537, 1, 2002; T. Padman- 
abhan, Phys. Rev. D, 66, 021301, 2002; T. Padmanab- 
han, T.R. Choudury, Phys. Rev. D, 66, 081301, 2002; J.S. 
Bagla, H.K. Jassal, T. Padmanabhan, Phys. Rev. D, 67, 
063504, 2003; E. Elizalde, S. Nojiri, S.D.Odintsov, Phys. 
Rev. D, 70, 040539, 2004 

[15] B.A. Bassett, M. Kunz, D. Parkinson, C. Ungarelli, Phys. 
Rev. D, 68, 043504 

[16] V.F. Cardone, A. Troisi, S. Capozziello, Phys. Rev. D, 

69, 083517, 2004 

[17] S. Capozziello, A. Melchiorri, A. Schirone, Phys. Rev. D, 

70, 101301, 2004 

[18] H. Sandvik et al., Phys. Rev. D, 69, 123524, 2004 
[19] A. Lue, R. Scoccimarro, G. Starkman, Phys. Rev. D, 69, 
044005, 2004 

[20] K. Freese, M. Lewis, Phys. Lett. B, 540, 1, 2002; K. 

Freese, Nucl. Phys. Suppl., 124, 50, 2003; Y. Wang, K. 

Freese, P. Gondolo, M. Lewis, ApJ, 594, 25, 2003 
[21] G.R. Dvali, G. Gabadadze, M. Porrati, Phys. Lett. B, 

485, 208, 2000; G.R. Dvali, G. Gabadadze, M. Kolanovic, 

F. Nitti, Phys. Rev. D, 64, 084004, 2001; G.R. Dvali, 

G. Gabadadze, M. Kolanovic, F. Nitti, Phys. Rev. D, 
64, 024031, 2002; A. Lue, R. Scoccimarro, G. Starkman, 
Phys. Rev. D, 69, 124015, 2004 

[22] R. Kerner, Gen. Rel. Grav. 14, 453, 1982 

[23] S. Capozziello, Int. J. Mod. Phys. D, 11, 483, 2002 

[24] S. Nojiri, S.D. Odintsov, Phys. Lett. B, 576, 5, 2003; S. 

Nojiri, S.D. Odintsov, Mod. Phys. Lett. A, 19, 627, 2003; 

S. Nojiri, S.D. Odintsov, Phys. Rev. D, 68, 12352, 2003; 

S.M. Carroll, V. Duvvuri, M. Trodden, M. Turner, Phys. 



15 



Rev. D, 70, 043528, 2004; S. Carloni, P.K.S. Dunsby, S. 
Capozziello, A. Troisi, astro - ph/0410046 
[25] S. Capozziello, S. Carloni, A. Troisi, astro- ph/0303041, 
2003 

[26] S. Nojiri, S.D. Odintsov, Gen. Rel. Grav., 36, 1765, 2004; 
X.H. Meng, P. Wang, Phys. Lett. B, 584, 1, 2004 

[27] I.L. Buchbinder, S.D. Odintsov, I.L. Shapiro, Effective 
Action in Quantum Gravity, IOP Publishing (1992) Bris- 
tol 

[28] C. Brans, R.H. Dicke, Phys. Rev., 124, 925 1961 

[29] N.D. Birrell, P.C.W. Davies, Quantum Fields in Curved 

Space, Cambridge Univ. Press, Cambridge (1982) 
[30] G. Vilkovisky, Class. Quant. Grav. 9, 895, 1992 
[31] G. Magnano, M. Ferraris, M. Francaviglia, Gen. Rel. 

Grav. 19, 465, 1987; M. Ferraris, M., Francaviglia, G. 

Magnano, Class. Quant. Grav. 7, 261, 1990 
[32] J. Barrow, A.C. Ottewill, J. Phys. A: Math. Gen., 16, 

2757, 1983 

[33] A. A. Starobinsky, Phys. Lett. B 91, 99 1980 
[34] D. La, P.J. Steinhardt, Phys. Rev. Lett. 62,376, 1989 
[35] P. Teyssandier, Ph. Tourrenc, J. Math. Phys. 24, 2793, 
1983 

[36] K. Maeda, Phys. Rev. D, 39, 3159, 1989 

[37] D. Wands, Class. Quant. Grav. 11, 269, 1994; S. 

Capozziello, R. de Ritis, A. A. Marino. Gen. Rel. Grav. 

30, 1247, 1998; S. Capozziello, G. Lambiase, Gen. Rel. 

Grav. 32, 295, 2000 
[38] S. Gottlober, H.-J. Schmidt, A. A. Starobinsky, Class. 

Quant. Grav. 7, 893, 1990 
[39] T.V. Ruzmaikina and A. A. Ruzmaikin, JETP, 30, 372, 

1970 

[40] L. Amendola, A. Battaglia-Mayer, S. Capozziello, S. 

Gottlober, V. Miiller, F. Occhionero, H.-J. Schmidt, 

Class. Quant. Grav. 10, L43, 1993 
[41] A. Battaglia-Mayer and H.-J. Schmidt, Class. Quant. 

Grav. 10, 2441, 1993 
[42] H.-J. Schmidt, Class. Quant. Grav. 7, 1023, 1990 
[43] L. Amendola, S. Capozziello, M. Litterio, F. Occhionero, 

Phys. Rev. D 45, 417, 1992; 
[44] CM. Will, Theory and Experiments in Gravitational 

Physics (1993) Cambridge Univ. Press, Cambridge 
[45] K. Stelle, Gen. Rel. Grav., 9, 353, 1978 
[46] R.H. Sanders, ARAA, 2 1, 1990 

[47] P.D. Mannheim, D. Kazanas, ApJ, 342, 635, 1989; O.V. 
Barabash, Yu. V. Shtanov, Phys. Rev. D, 60, 064008, 
1999 

[48] J.D. Anderson et al., Phys. Rev. Lett. 81, 2858, 1998; 

J.D. Anderson, et al., Phys. Rev. D, 65, 082004, 2002 
[49] I. Quant, H.-J. Schmidt, Astron. Nachr. 312, 97, 1991 
[50] P. Schneider, J. Ehlers,E.E. Falco, Gravitational Lenses, 

Springer- Verlag (1992) Berlin. 
[51] L.M. Krauss, M. White, ApJ. 397, 357, 1992 



[52] L. Nottale in Dark Matter (Moriond Astrophysics Meet- 
ings), J. Andouze and J. Tran Thanh Van eds. (1988) 
Frontieres, Gif-sur-Yvette. 

[53] M. Ferraris, M. Francaviglia, I. Volovich, Class. Quant. 
Grav., 11, 1505, 1994 

[54] G. Allemandi, A. Borowiec, M. Francaviglia, Phys. Rev. 
D, 70, 043524, 2004; G. Allemandi, A. Borowiec, M. 
Francaviglia, Phys. Rev. D, 70, 103503, 2004; 

[55] G. Allemandi, M. Capone, S. Capozziello, M. Francav- 
iglia, hep-th/0409198, 2004 

[56] D.N. Vollick, Phys. Rev. 68, 063510, 2003; X.H. Meng, P. 
Wang, Class. Quant. Grav., 20, 4949, 2003; E.E. Flana- 
gan, Phys. Rev. Lett. 92, 071101, 2004; E.E. Flanagan, 
Class. Quant. Grav., 21, 417, 2004; X.H. Meng, P. Wang, 
Class. Quant. Grav., 21, 951, 2004 G.M. Kremer, D.S.M. 
Alves, Phys. Rev. D, 70, 023503, 2004 

[57] S. Capozziello, V.F. Cardone, S. Carloni, A. Trosi, Int. 
J. Mod. Phys. D 12, 1969, 2003 

[58] J. P. Duirisseau, R. Kerner, P. Eysseric, Gen. Rel. Grav., 
15, 797, 1983; J. P. Duirisseau, R. Kerner, Class. Quant. 
Grav. 3, 817, 1986 

[59] M. Visser, Class. Quant. Grav., 21, 2603, 2004 

[60] P.S. Corasaniti, M. Kunz, D. Parkinson, E.J. Copeland 
and B.A. Bassett, astro - ph/0406608, 2004 

[61] V. Sahni, T.D. Saini, A. A. Starobinsky, U. Alam, JETP 
Lett., 77, 201, 2003 

[62] R. Caldwell, M. Kamionkowski, JCAP, 09, 9, 2004; M.V. 
John, ApJ, 614, 1, 2004 

[63] A. A. Starobinksy, JETP Lett., 68, 757, 1998; T. Naka- 
mura, T. Chiba, MNRAS, 306, 696, 1999 

[64] D. Huterer, M.S. Turner, Phys. Rev. D, 081301, 1999 

[65] R. Dick, Gen. Rel. Grav., 36, 217, 2004; A.E. Dominguez, 
D.E. Barraco, Phys. Rev. D, 70, 043505, 2004 

[66] T.D. Saini, S. Raychauduri, V. Sahni, A. A. Starobinski, 
Phys. Rev. Lett., 85, 1162, 2000 

[67] U. Alam, V. Sahni, T.D. Saini, A. A. Starobinsky, MN- 
RAS, 344, 1057, 2003; U. Alam, V. Sahni, T.D. Saini, 
A. A. Starobinsky, MNRAS, 354, 275, 2004; U. Alam, V. 
Sahni, A. A. Starobinsky, JCAP, 0406, 008, 2004; J. Jon- 
sson, A. Goobar, R. Amanullah, L. Bergstrom, JCAP, 
0409, 007, 2004; U. Alam, V. Sahni, T.D. Saini, A. A. 
Starobinsky, astro- ph/0406672, 2004 

[68] S. Hannestad, E. Mortsell, Phys. Rev. D, 66, 0635088, 
2002; A. Melchiorri, L. Mersini, C.J. Odman, M. Trod- 
den, Phys. Rev. D, 68, 043509, 2003; S. Hannestad, E. 
Mortsell, JCAP, 0409, 001, 2004 

[69] G. Allemandi, S. Capozziello, V.F. Cardone, M. Francav- 
iglia, A. Troisi, in preparation. 

[70] A.D. Dolgov, M. Kawasaki, Phys. Lett. B, 573, 1, 2003 

[71] S. Nojiri, astro -ph/0407099, 2004 



